<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.5"/>
<title>mechanoChemIGA: Example 2 : Cahn-Hilliard (one species)</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
  $(document).ready(function() { searchBox.OnSelectItem(0); });
</script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  <td style="padding-left: 0.5em;">
   <div id="projectname">mechanoChemIGA
   &#160;<span id="projectnumber">0.2.2</span>
   </div>
  </td>
 </tr>
 </tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.5 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
  <div id="navrow1" class="tabs">
    <ul class="tablist">
      <li><a href="index.html"><span>Main&#160;Page</span></a></li>
      <li><a href="codestructure.html"><span>Code&#160;structure</span></a></li>
      <li><a href="examples.html"><span>Examples</span></a></li>
      <li><a href="modules.html"><span>Functions</span></a></li>
      <li><a href="annotated.html"><span>Classes</span></a></li>
      <li><a href="files.html"><span>Files</span></a></li>
      <li><a href="https://github.com/mechanoChem/mechanoChem"><span>GitHub&#160;code&#160;repo</span></a></li>
      <li>
        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.png"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
          </span>
        </div>
      </li>
    </ul>
  </div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Variables</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark">&#160;</span>Typedefs</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(6)"><span class="SelectionMark">&#160;</span>Enumerations</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(7)"><span class="SelectionMark">&#160;</span>Enumerator</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(8)"><span class="SelectionMark">&#160;</span>Macros</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(9)"><span class="SelectionMark">&#160;</span>Groups</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(10)"><span class="SelectionMark">&#160;</span>Pages</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

</div><!-- top -->
<div class="header">
  <div class="headertitle">
<div class="title">Example 2 : Cahn-Hilliard (one species) </div>  </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><p>This example implements the Cahn-Hilliard equation for phase-field modeling of a single species, as described by the following weak form of the PDE. The scalar field is composition, <img class="formulaInl" alt="$c$" src="form_5.png"/>. Note the application of the higher-order Dirichlet boundary condition <img class="formulaInl" alt="$\nabla c\cdot\boldsymbol{n}=0$" src="form_7.png"/> using Nitsche's method.</p>
<p>Cahn-Hilliard:</p>
<p class="formulaDsp">
<img class="formulaDsp" alt="\begin{eqnarray*} 0 &amp;=&amp; \int_\Omega \left(w_1\frac{c - c_{prev}}{\mathrm{d}t} + M\left(\nabla w_1\cdot(f_{,cc}\nabla c) + \kappa_1\nabla^2 w_1\nabla^2 c\right)\right) dV\\ &amp;\phantom{=}&amp; - \int_{\partial\Omega} \left(w_1j_n + M\kappa_1\left(\nabla^2c(\nabla w_1\cdot\boldsymbol{n}) + \nabla^2w_1(\nabla c\cdot\boldsymbol{n})\right) - \tau(\nabla w_1\cdot\boldsymbol{n})(\nabla c\cdot\boldsymbol{n})\right) dS \end{eqnarray*}" src="form_24.png"/>
</p>
<p>Free energy density:</p>
<p class="formulaDsp">
<img class="formulaDsp" alt="\begin{eqnarray*} f(c) = \alpha(c - c_a)^2(c - c_b)^2 \end{eqnarray*}" src="form_25.png"/>
</p>
<p>The current settings prescribe random initial conditions and zero-flux boundary conditions. With these settings, the following evolution of the concentration is obtained:</p>
 <style>div.image img[src="example2.png"]{width:20cm;}</style>  <div class="image">
<img src="example2.png" alt="example2.png"/>
</div>
<h1>Implementation: Level 1 users </h1>
<p>To implement this model, we will specify the following through defining user functions: <br/>
</p>
<ul>
<li>Initial conditions <br/>
</li>
<li>Constitutive model (via free energy density functions) <br/>
</li>
<li>Parameter values <br/>
</li>
<li>Weak form of the PDEs <br/>
</li>
</ul>
<p>First, we include the header file declaring the required user functions. These functions will be defined in this file.</p>
<p><div class="fragment"><div class="line"><span class="preprocessor">#include &quot;<a class="code" href="user_functions_8h.html">userFunctions.h</a>&quot;</span></div>
</div><!-- fragment --></p>
<p>Now, we first define any optional user functions. Optional user functions have a default definition that can be redefined by the user using a function pointer. This will be done in the <code>defineParameters</code> function. The available list of optional user functions includes: <code>boundaryConditions</code>, <code>scalarInitialConditions</code>, <code>vectorInitialConditions</code>, <code>loadStep</code>, <code>adaptiveTimeStep</code>, and <code>projectFields</code>. In this example, we redefine only the <code>scalarInitialConditions</code> function, while using the default functions for the others.</p>
<p><b> The <code>scalarInitialConditions</code> function </b></p>
<p>We initialized the composition field to be random about 0.5.</p>
<p><div class="fragment"><div class="line"><span class="keyword">template</span>&lt;<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> dim&gt;</div>
<div class="line"><span class="keywordtype">double</span> <a class="code" href="_cahn_hilliard___allen_cahn_22_d_2user_functions_8cc.html#a360bfcfddf6eb0952d9904c883f3f741">userScalarInitialConditions</a>(<span class="keyword">const</span> <a class="code" href="class_tensor.html">Tensor&lt;1,dim,double&gt;</a> &amp;x, <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> scalar_i, <a class="code" href="struct_app_ctx.html">AppCtx&lt;dim&gt;</a> &amp;user)</div>
<div class="line">{</div>
<div class="line"></div>
<div class="line">  <span class="keywordflow">return</span> 0.5 + 0.3*(0.5 - (double)(rand() % 100 )/100.0); <span class="comment">//Random about 0.5</span></div>
<div class="line"></div>
<div class="line">} <span class="comment">//end scalarInitialConditions</span></div>
</div><!-- fragment --></p>
<p><b> Free energy density derivative functions </b></p>
<p>This phase-field implementation requires the second derivative of the chemical free energy density function <img class="formulaInl" alt="$f(c,\eta) = \alpha(c - c_a)^2(c - c_b)^2$" src="form_26.png"/>. We define the function computing <img class="formulaInl" alt="$\partial^2 f/\partial c^2$" src="form_43.png"/> here. Note that this free energy derivative function is used only in this file. It is not a member of any class, nor will we use it to set any function pointers.</p>
<p><div class="fragment"><div class="line"><span class="keyword">template</span>&lt;<span class="keyword">typename</span> T&gt;</div>
<div class="line">T <a class="code" href="_cahn_hilliard___allen_cahn_22_d_2user_functions_8cc.html#acf9953d272a25dde719af6ce9d2763cb">F_cc</a>(T c, <span class="keywordtype">double</span> alpha, <span class="keywordtype">double</span> ca, <span class="keywordtype">double</span> cb){</div>
<div class="line">  </div>
<div class="line">  <span class="comment">//Second derivative of the free energy density: f(c) = alpha*(c-ca)^2*(c-cb)^2</span></div>
<div class="line">  <span class="keywordflow">return</span> 2*alpha*(std::pow(2.*c-ca-cb,2) + 2.*(c-ca)*(c-cb));</div>
<div class="line">  </div>
<div class="line">} <span class="comment">//end F_cc</span></div>
</div><!-- fragment --></p>
<p><b> The <code>defineParameters</code> function </b></p>
<p>The user is required to define the <code>defineParameters</code> and <code>residual</code> functions. The <code>defineParameters</code> defines variables and functions in the <code>AppCtx</code> object. The <code>AppCtx</code> object is defined in the appCtx.h file. This function is used to define any values in <code>user</code> that will be needed in the problem. It is also used to set any function pointers for user functions that we have redefined.</p>
<p>Many of these values can be overwritten by the parameters.prm file, which we will look at later.</p>
<p><div class="fragment"><div class="line"><span class="keyword">template</span>&lt;<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> dim&gt;</div>
<div class="line"><span class="keywordtype">void</span> <a class="code" href="group__user_functions.html#gadbccf6631ad847d5a681a548f921ef29">defineParameters</a>(<a class="code" href="struct_app_ctx.html">AppCtx&lt;dim&gt;</a>&amp; user){</div>
</div><!-- fragment --></p>
<p>Here, we define the mesh by setting the number of elements in each direction, e.g. a 100x100 element mesh.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a55ccabd543df9a0223cd34dbd64c987d">N</a>[0] = 100;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a55ccabd543df9a0223cd34dbd64c987d">N</a>[1] = 100;</div>
</div><!-- fragment --></p>
<p>We also define the dimensions of the domain, e.g. a unit square.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a789652912f4d6df6c0836aa22ae93de0">L</a>[0] = 1.;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a789652912f4d6df6c0836aa22ae93de0">L</a>[1] = 1.;</div>
</div><!-- fragment --></p>
<p>We can define a periodic (or partially periodic) domain. The default is no periodicity in all directions. Here, we override the default and define periodicity in the x direction.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#ac77e8a2af239ac3f36a49dfc69ad5c69">periodic</a>[0] = PETSC_TRUE;</div>
</div><!-- fragment --></p>
<p>We can define additional material parameters that are not explicity listed in the <code>user</code> structure by defining elements of the <code>matParam</code> C++ map, which maps <code>std::string</code> to <code>double</code>. These values can also be overwritten in the parameters file.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;inFlux&quot;</span>] = 0; <span class="comment">//Infux through top</span></div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;mobility&quot;</span>] = .1; <span class="comment">//Mobility</span></div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;kappa&quot;</span>] = .0001; <span class="comment">//Gradient energy parameter</span></div>
<div class="line"></div>
<div class="line">  <span class="comment">//Define some free energy parameters</span></div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;alpha&quot;</span>] = 0.25; <span class="comment">//Free energy coefficient</span></div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;c_a&quot;</span>] = 0.2; <span class="comment">//Composition of phase a</span></div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;c_b&quot;</span>] = 0.9; <span class="comment">//Composition of phase b</span></div>
</div><!-- fragment --></p>
<p>We define the initial time step and total simulation time. We also have the options to use restart files, in which case we would set the iteration index and time at which to start. We leave these values at zero to begin a new simulation. We also have the option to output results at regular intervals (e.g. every 5 time steps).</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a4c155b92216444548c4457f18e050630">dtVal</a> = .1;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a17b4070be131d7cee7880bbc695f9169">totalTime</a> = 20;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a28a6be93b52da95fb5e891f10a9a5d87">RESTART_IT</a> = 0;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a6d5465cb515bc053db8785edd2e1b21b">RESTART_TIME</a> = 0.;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#a7c81203929d679b6ee72dec416c9bea9">skipOutput</a> = 2;</div>
</div><!-- fragment --></p>
<p>We specify the number of vector and scalar solution and projection fields by adding the name of each field to their respective vector. Here, we have one scalar solution field (the composition). We do not use any vector solution fields or projection fields in this example.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a6794fa3a8512e9239c9f2e4c7bf18ec3">scalarSolnFields</a>.push_back(<span class="stringliteral">&quot;c&quot;</span>);</div>
</div><!-- fragment --></p>
<p>We can specify the polynomial order of the basis splines, as well as the global continuity. Note that the global continuity must be less than the polynomial order. Here, we use quadratic basis functions with C-1 global continuity.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a23cded74ca3d8ec2f99d69e41d8539ca">polyOrder</a> = 2;</div>
<div class="line">  user.<a class="code" href="struct_app_ctx.html#ae091f872a8ec5d2a5492586ae8fcbcbb">globalContinuity</a> = 1;</div>
</div><!-- fragment --></p>
<p>Finally, we redirect the desired user function pointers to the <code>scalarInitialConditions</code> function that we defined above. This completes the <code>defineParameters</code> function.</p>
<p><div class="fragment"><div class="line">  user.<a class="code" href="struct_app_ctx.html#a4a4f8fca81419b6f43a612a72ed0a989">scalarInitialConditions</a> = <a class="code" href="_cahn_hilliard___allen_cahn_22_d_2user_functions_8cc.html#a360bfcfddf6eb0952d9904c883f3f741">userScalarInitialConditions</a>;</div>
<div class="line"></div>
<div class="line">} <span class="comment">//end defineParameters</span></div>
</div><!-- fragment --></p>
<p><b> The <code>residual</code> function </b></p>
<p>The residual function defines the residual that is to be driven to zero. This is the central function of the code. It is set up to follow the analytical weak form of the PDE. It has a number of arguments that give problem information at the current quadrature point.</p>
<p><div class="fragment"><div class="line"><span class="keyword">template</span>&lt;<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> dim, <span class="keyword">typename</span> T&gt;</div>
<div class="line"><span class="keywordtype">void</span> <a class="code" href="group__user_functions.html#gab9195b3f02c923dafb2c742df293db7d">residual</a>(<span class="keywordtype">bool</span> dV,</div>
<div class="line">          <span class="keywordtype">bool</span> dS,</div>
<div class="line">          <span class="keyword">const</span> <a class="code" href="class_tensor.html">Tensor&lt;1,dim,double&gt;</a> &amp;x,</div>
<div class="line">          <span class="keyword">const</span> <a class="code" href="class_tensor.html">Tensor&lt;1,dim,double&gt;</a> &amp;normal,</div>
<div class="line">          <span class="keyword">const</span> <a class="code" href="classsolution_scalars.html">solutionScalars&lt;dim,T&gt;</a> &amp;c,</div>
<div class="line">          <span class="keyword">const</span> <a class="code" href="classsolution_vectors.html">solutionVectors&lt;dim,T&gt;</a> &amp;u,</div>
<div class="line">          <span class="keyword">const</span> <a class="code" href="classtest_function_scalars.html">testFunctionScalars&lt;dim,T&gt;</a> &amp;w1,</div>
<div class="line">          <span class="keyword">const</span> <a class="code" href="classtest_function_vectors.html">testFunctionVectors&lt;dim,T&gt;</a> &amp;w2,</div>
<div class="line">          <a class="code" href="struct_app_ctx.html">AppCtx&lt;dim&gt;</a> &amp;user,</div>
<div class="line">          Sacado::Fad::SimpleFad&lt;T&gt; &amp;r){</div>
</div><!-- fragment --></p>
<p><code>dV</code> is a boolean, "true" if <code>residual</code> is being called for the volume integral and "false" if <code>residual</code> is being called for the surface integral.<br/>
<code>dS</code> is a boolean, "false" if <code>residual</code> is being called for the volume integral and "true" if <code>residual</code> is being called for the surface integral.<br/>
<code>x</code> gives the coordinates of the quadrature point.<br/>
<code>normal</code> gives the unit normal for a surface quadrature point.<br/>
<code>c</code> gives the information (values, gradients, etc.) for the scalar solution fields at the current quadrature point (see documentation for solutionScalars class).<br/>
<code>u</code> gives the information (values, gradients, etc.) for the vector solution fields at the current quadrature point (see documentation for solutionVectors class).<br/>
<code>w1</code> gives the information for the scalar test functions.<br/>
<code>w2</code> gives the information for the vector test functions.<br/>
<code>user</code> is a structure available for parameters related to the initial boundary value problem (e.g. elasticity tensor).<br/>
<code>r</code> stores the scalar value of the residual for the weak form of the PDE which is then used by the core assembly functions.</p>
<p>The following functions are available for the solution objects <code>c</code> and <code>u</code>, where the argument is the field index, i.</p>
<p><code>c.val(i)</code> - Value of scalar field i, scalar <br/>
<code>c.grad(i)</code> - Gradient of scalar field i, 1st order tensor <br/>
<code>c.hess(i)</code> - Hessian of scalar field i, 2nd order tensor <br/>
<code>c.laplacian(i)</code> - Laplacian of scalar field i, scalar <br/>
<code>c.valP(i)</code> - Value of scalar field i at previous time step, scalar <br/>
<code>c.gradP(i)</code> - Gradient of scalar field i at previous time step, 1st order tensor <br/>
<code>c.hessP(i)</code> - Hessian of scalar field i at previous time step, 2nd order tensor <br/>
<code>c.laplacianP(i)</code> - Laplacian of scalar field i at previous time step, scalar</p>
<p><code>u.val(i)</code> - Value of vector field i, 1st order tensor <br/>
<code>u.grad(i)</code> - Gradient of vector field i, 2nd order tensor <br/>
<code>u.hess(i)</code> - Hessian of vector field i, 3rd order tensor <br/>
<code>u.valP(i)</code> - Value of vector field i at previous time step, 1st order tensor <br/>
<code>u.gradP(i)</code> - Gradient of vector field i at previous time step, 2nd order tensor <br/>
<code>u.hessP(i)</code> - Hessian of vector field i at previous time step, 3rd order tensor</p>
<p>Similar functions are available for the test functions. Also, the following tensor operations are useful:</p>
<p>Tensor operations: <br/>
<code>operator+</code> - tensor addition <br/>
<code>operator-</code> - tensor subraction <br/>
<code>operator*</code> - single contraction between tensors or scalar multiplication <br/>
<code>double_contract</code> - double contraction of two 2nd order tensors, or a 4th order tensor and a 2nd order tensor. <br/>
<code>trans( )</code> - transpose 2nd order tensor <br/>
<code>trace( )</code> - trace of 2nd order tensor <br/>
<code>det( )</code> - determinant of 2nd order tensor <br/>
<code>inv( )</code> - inverse of 2nd order tensor <br/>
</p>
<p>The example code here implements the weak form for the Cahn-Hilliard equation, as shown above.</p>
<p>First, we set the values for necessary parameters, using some predefined material parameters.</p>
<p><div class="fragment"><div class="line">  <span class="keywordtype">double</span> dt = user.<a class="code" href="struct_app_ctx.html#a102ae36ef6e9df9f12c7e988659b5ff0">dt</a>;</div>
<div class="line">  <span class="keywordtype">double</span> jn = user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;inFlux&quot;</span>]*(x[1]==user.<a class="code" href="struct_app_ctx.html#a789652912f4d6df6c0836aa22ae93de0">L</a>[1]); <span class="comment">//Influx through the top</span></div>
<div class="line">  <span class="keywordtype">double</span> M = user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;mobility&quot;</span>]; <span class="comment">//Mobility</span></div>
<div class="line">  <span class="keywordtype">double</span> kappa = user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;kappa&quot;</span>];</div>
<div class="line">  <span class="keywordtype">double</span> tau = 0.1*(user.<a class="code" href="struct_app_ctx.html#a55ccabd543df9a0223cd34dbd64c987d">N</a>[0]/user.<a class="code" href="struct_app_ctx.html#a789652912f4d6df6c0836aa22ae93de0">L</a>[0]);</div>
</div><!-- fragment --></p>
<p>Next, we get the values for the free energy derivative <img class="formulaInl" alt="$f_{,cc}$" src="form_16.png"/> based on the current quadrature point.</p>
<p><div class="fragment"><div class="line">  T f_cc;</div>
<div class="line">  f_cc = <a class="code" href="_cahn_hilliard___allen_cahn_22_d_2user_functions_8cc.html#acf9953d272a25dde719af6ce9d2763cb">F_cc</a>(c.<a class="code" href="classsolution_scalars.html#aad0d4f18ccffe03f32a1ac985eb3b3d8">val</a>(0),user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;alpha&quot;</span>],user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;c_a&quot;</span>],user.<a class="code" href="struct_app_ctx.html#a36448fa26553468f32e3cab046566bea">matParam</a>[<span class="stringliteral">&quot;c_b&quot;</span>]);</div>
</div><!-- fragment --></p>
<p>Now, we compute the residual in a manner very similar to the analytical form:</p>
<p class="formulaDsp">
<img class="formulaDsp" alt="\begin{eqnarray*} 0 &amp;=&amp; \int_\Omega \left(w_1\frac{c - c_{prev}}{\mathrm{d}t} + M\left(\nabla w_1\cdot(f_{,cc}\nabla c) + \kappa_1\nabla^2 w_1\nabla^2 c\right)\right) dV\\ &amp;\phantom{=}&amp; - \int_{\partial\Omega} \left(w_1j_n + M\kappa_1\left(\nabla^2c(\nabla w_1\cdot\boldsymbol{n}) + \nabla^2w_1(\nabla c\cdot\boldsymbol{n})\right) - \tau(\nabla w_1\cdot\boldsymbol{n})(\nabla c\cdot\boldsymbol{n})\right) dS \end{eqnarray*}" src="form_24.png"/>
</p>
<p><div class="fragment"><div class="line">  r = ( w1.<a class="code" href="classtest_function_scalars.html#ad76f9644680c80e6602ccf460600b3a3">val</a>(0)*(c.<a class="code" href="classsolution_scalars.html#aad0d4f18ccffe03f32a1ac985eb3b3d8">val</a>(0) - c.<a class="code" href="classsolution_scalars.html#aa4c4b85dc344e177b1cbeca4bd59f0e8">valP</a>(0))/dt )*dV;  </div>
<div class="line">  r += M*w1.<a class="code" href="classtest_function_scalars.html#a40dbbf7383a6a9471d897f455405315d">grad</a>(0)*f_cc*c.<a class="code" href="classsolution_scalars.html#a94a7f34ce00b9f6c3584756b185b320b">grad</a>(0)*dV;</div>
<div class="line">  r += M*kappa*w1.<a class="code" href="classtest_function_scalars.html#a97d1e6abaec05fc0971181ad38af50fd">laplacian</a>(0)*c.<a class="code" href="classsolution_scalars.html#aaede18812f64934bc57ba259320d1240">laplacian</a>(0)*dV;</div>
<div class="line">  </div>
<div class="line">  r += -w1.<a class="code" href="classtest_function_scalars.html#ad76f9644680c80e6602ccf460600b3a3">val</a>(0)*jn*dS; <span class="comment">//Boundary flux</span></div>
<div class="line">  r += -M*kappa*( c.<a class="code" href="classsolution_scalars.html#aaede18812f64934bc57ba259320d1240">laplacian</a>(0)*(w1.<a class="code" href="classtest_function_scalars.html#a40dbbf7383a6a9471d897f455405315d">grad</a>(0)*normal) + w1.<a class="code" href="classtest_function_scalars.html#a97d1e6abaec05fc0971181ad38af50fd">laplacian</a>(0)*(c.<a class="code" href="classsolution_scalars.html#a94a7f34ce00b9f6c3584756b185b320b">grad</a>(0)*normal) )*dS;</div>
<div class="line">  r += tau*(w1.<a class="code" href="classtest_function_scalars.html#a40dbbf7383a6a9471d897f455405315d">grad</a>(0)*normal)*(c.<a class="code" href="classsolution_scalars.html#a94a7f34ce00b9f6c3584756b185b320b">grad</a>(0)*normal)*dS;</div>
<div class="line">  </div>
<div class="line">} <span class="comment">//end residual</span></div>
</div><!-- fragment --></p>
<p>Finally, we include a file that instatiates the template functions <code>defineParameters</code> and <code>residual</code>. This bit of code will generally be the same for any problem (unless you decide to use a different automatic differentation library); the user does not need to modify it.</p>
<p><div class="fragment"><div class="line"><span class="preprocessor">#include &quot;<a class="code" href="user_functions_instantiation_8h.html">userFunctionsInstantiation.h</a>&quot;</span></div>
</div><!-- fragment --></p>
<p>The complete implementation can be found at <a href="https://github.com/mechanoChem/mechanoChemIGA/blob/master/initBounValProbs/CahnHilliard_oneSpecies/2D/userFunctions.cc">Github</a>.</p>
<h1>Parameters file: Interface for level 2 users </h1>
<p>Now let's look at the parameters file, <code>parameters.prm</code>. The advantages of the parameters file are that these values can be changed without recompiling the code and it can provide a clean interface to the code. <p>The parameters defined in the parameters file overwrite any previous values defined in the <code>defineParameters</code> function. Anything following the pound sign (#) is a comment. A parameter is defined using the syntax:</p>
<p><code>set</code> <code>parameterName</code> <code>=</code> <code>parameterValue</code> </p>
<p>There is a set list of variables that can be read from the parameters file. Anything else will be added to the <code>matParam</code> structure with a double number type. Tensor objects can follow the format: 1 x 1 or [1,1] or (1,1), where the number of components must equal the spatial dimension of the problem.</p>
<p>In this example file, we begin by specifying the spatial dimension, the geometry dimensions, and the mesh size:</p>
<p><div class="fragment"><div class="line">set dim = 2  # spatial dimension</div>
<div class="line">set L = 1 x 1  # dimensions of geometry</div>
<div class="line">set N = 100 x 100 # mesh discretization</div>
</div><!-- fragment --></p>
<p>Next, we define some parameters that are specific to this problem, so they become elements of <code>matParam</code> (see the <code>residual</code> and  functions above).</p>
<p><div class="fragment"><div class="line"><span class="preprocessor"># Free energy parameters, f(c) = alpha*(c-c_a)^2*(c-c_b)^2</span></div>
<div class="line"><span class="preprocessor"></span>set alpha = 0.25  # free energy coefficient</div>
<div class="line">set c_a = 0.2  # composition of phase a</div>
<div class="line">set c_b = 0.7  # composition of phase b</div>
<div class="line"></div>
<div class="line"><span class="preprocessor"># Material parameters</span></div>
<div class="line"><span class="preprocessor"></span>set inFlux = 0.  # flux through the top boundary</div>
<div class="line">set mobility = 0.1</div>
<div class="line">set kappa = 0.0001  # gradient energy parameter</div>
</div><!-- fragment --></p>
<p>We then define time stepping, restart information, output frequency, and spline parameters.</p>
<p><div class="fragment"><div class="line"><span class="preprocessor"># Time stepping</span></div>
<div class="line"><span class="preprocessor"></span>set dtVal = 1.  # initial time step size</div>
<div class="line">set totalTime = 100.</div>
<div class="line"></div>
<div class="line"># Restart</div>
<div class="line">set RESTART_TIME = 0.</div>
<div class="line">set RESTART_IT = 0</div>
<div class="line"></div>
<div class="line"><span class="preprocessor"># Output</span></div>
<div class="line"><span class="preprocessor"></span>set skipOutput = 1</div>
<div class="line"></div>
<div class="line"><span class="preprocessor"># Splines</span></div>
<div class="line"><span class="preprocessor"></span>set polyOrder = 2</div>
<div class="line">set globalContinuity = 1</div>
</div><!-- fragment --></p>
<p>Note that we don't need to include all (or even any) of these parameters in this file. We defined default values previously.</p>
<p>The complete parameters file can be found at <a href="https://github.com/mechanoChem/mechanoChemIGA/blob/master/initBounValProbs/CahnHilliard_oneSpecies/2D/parameters.prm">Github</a>. </p>
</div></div><!-- contents -->
<!-- start footer part -->
<hr class="footer"/><address class="footer"><small>
Generated by &#160;<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.8.5
</small></address>
</body>
</html>
